#include "initialSol.h"
using namespace std;

//------------------------------------------------------------------------------
InitialSolution::InitialSolution(Mesh* mymesh,FsContainer* myfs,SaveOutputGmsh* output,ClassParameter& param){
  
  mesh        = mymesh;
  fs          = myfs; 
  initSolName = param.getToken("InitialSol");
  
  if(initSolName=="Smooth")             initTypeSol=0;
  if(initSolName=="ThreeBodyRotation")  initTypeSol=1;
  if(initSolName=="Burger")             initTypeSol=2;
  if(initSolName=="KPP")                initTypeSol=3;
  if(initSolName=="SineWave")           initTypeSol=4;

  // --- initialize Uinit container ---
  
  int eT=0;
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt,eT++){
    Uinit.push_back(DGMatrix((*fsIt)->getNbENP(),mymesh->getNbElmt(eT)));
  }

  // --- compute initial sol ---

  computeInitialSolution();
  
  // --- save output file ---

  output->saveOutput(&Uinit,0,0,1,"solution");
  

}
//------------------------------------------------------------------------------
double InitialSolution::analyticalSol(double x,double y){
  
  switch (initTypeSol) {
    case Smooth:            return SmoothSolution(x,y);break;
    case ThreeBodyRotation: return ThreeBodyRotationSolution(x,y);break;
    case Burger:            return BurgerSolution(x,y);break;
    case KPP:               return KPPSolution(x,y);break;
    case SineWave:          return sin(2*M_PI*x/1.);break;
    default:                return 1.;break;
  }

}
//------------------------------------------------------------------------------
void InitialSolution::computeInitialSolution(){
  
  int eT=0;
  for(FsContainer::const_iterator fsIt = fs->begin();fsIt!=fs->end();++fsIt,eT++){
    
    const DGMatrix* ERefNodesPos = (*fsIt)->getERefNodesPos();
    
    FunctionSpace* mapping = FunctionSpaceFactory::CreateNewFunctionSpace((FctType)mesh->getEGmshType(eT),1);
    
    int nbENP = (*fsIt)->getNbENP();
    int nbF   = mesh->getNbElmtFaces(eT);

    #pragma omp parallel for
    for(int eID=0;eID<mesh->getNbElmt(eT);eID++){
      DGMatrix NPcoordXY(nbENP,2);
      double detJ = mesh->getDetJ(0,eID);
      for(int n=0;n<nbENP;n++){
        for(int j=0;j<nbF;j++){
          double phi_j    = mapping->getPhiE(j,(*ERefNodesPos)(n,0),(*ERefNodesPos)(n,1)); 
          NPcoordXY(n,0) += mesh->getNode(mesh->getElmtNodesTag(eT,eID,j),0)*phi_j;
          NPcoordXY(n,1) += mesh->getNode(mesh->getElmtNodesTag(eT,eID,j),1)*phi_j;
        }
      }  
      for(int i=0;i<nbENP;i++) Uinit[eT](i,eID) = analyticalSol(NPcoordXY(i,0),NPcoordXY(i,1));
    }
    
    delete mapping;
    
  }
  
}
//==============================================================================
//------------------------------------------------------------------------------
double InitialSolution::SmoothSolution(double x,double y){
  //dom = [-1,+1]^2
  
  double x_0[2] = {0.3,0.};
  double r_0    = 0.25;
  double dist   = pow(x-x_0[0],2)+pow(y-x_0[1],2);
  
  return  0.5*(1-tanh(dist/(r_0*r_0)-1));
  
}
//------------------------------------------------------------------------------
double InitialSolution::ThreeBodyRotationSolution(double x,double y){
  //dom = [-1,+1]^2
  
  double x_d[2] = {0.,0.5};
  double x_c[2] = {0.,-0.5};
  double x_h[2] = {-0.5,0.};
  double r_0    = 0.3;
  double norm_d = sqrt(pow(x-x_d[0],2)+pow(y-x_d[1],2));
  double norm_c = sqrt(pow(x-x_c[0],2)+pow(y-x_c[1],2));
  double norm_h = sqrt(pow(x-x_h[0],2)+pow(y-x_h[1],2));
  
  if(norm_d<=r_0 && (fabs(x)>=0.05 || y>=0.7)){
    return 1.;
  }else if(norm_c <= r_0){
    return 1-norm_c/r_0;
  }else if(norm_h<=r_0){
    return 0.5*(1+cos(M_PI*min(norm_h/r_0,1.)));
  }
  
  return 0.;
  
}
//------------------------------------------------------------------------------
double InitialSolution::BurgerSolution(double x,double y){
  //dom = [0,+1]^2
  
  if(x<=0.5){
    if(y>=0.5){
      return -0.2;
    }else{
      return  0.5;
    }
  }else{
    if(y>=0.5){
      return -1.;
    }else{
      return  0.8;
    }
  }
  
  return 0;
}
//------------------------------------------------------------------------------
double InitialSolution::KPPSolution(double x,double y){
  //dom = [-2,2]x[-2.5,1.5]
  
  if(sqrt(x*x+y*y)<=1){
    return 14*M_PI/4.;
  }else{
    return M_PI/4.;
  }
  
}